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ON  THE  SIMULATION  OF  HARMONICALLY  RELATED  SIGNALS 


INTRODUCTION 

In  a  number  of  applications,  narrowband  signal  components  occur  that  are  harmonically  related. 
To  understand  harmonically  related  signals,  it  is  sufficient  to  recognize  that  over  any  observation  (or 
analysis)  interval,  a  narrowband  signal  u(t)  may  be  represented  as  A(t)  exp  where  A{t)  >  0 

is  the  amplitude  function  and  \p(t)  is  the  phase  function  of  the  relevant  signal.  Any  two  narrowband 
signals  up(t )  and  uq (t )  are  considered  to  be  harmonically  related  when  the  phase  functions  of  the  two 
signals  are  linearly  related.  That  is,  when 

\p  (t)  =  qip(t)  +  a  constant  ( 1  a) 

and  the  harmonic  ratio  q  is  any  real  number  (not  necessarily  an  integer).  The  instantaneous  frequen¬ 
cies  of  the  two  functions  are  therefore  related  as 

/<?(')  =  Mf)/ 2lr  =  fq  +  Vq (f) 


=  q$r(t)/2ic  =  q[fp  +  vp(t)]  =  qfp(t ),  (lb) 

where  the  "dot”  over  the  variable  implies  the  time  derivative.  Here  fp(t),  fq(t)  are  the  instantane¬ 
ous  frequencies  of  the  two  signals,  fp,  fq  are  the  mean  values  of  the  instantaneous  frequencies,  and 
vp(t),  vq(t)  are  the  zero-mean  fluctuating  components.  For  convenience,  it  is  assumed  that 
fpU)  <  fgU)  and  1  <  q. 

This  report  formulates  discrete  mathematical  functions  that  simulate  the  harmonic  relations  given 
in  Eq.  (1)  for  implementation  on  a  digital  computer.  The  fluctuating  components  are  random  but 
bounded  in  their  excursion,  and  the  fluctuation  dynamics  (rate  of  change)  are  controlled  to  emulate 
signals  that  occur  in  practice. 
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FORMULATION  OF  THE  PHASE  FUNCTIONS 

To  formulate  the  discrete  phase  functions  of  harmonically  related  signals  sampled  at  a  rate  of  /, 
Hz,  let 
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This  expression  may  be  written  in  recurrence  form  as 


Z'n  =  z'n_,  + 

J  S 

where  Z'  _x  =0.  Since  phase  need  be  determined  only  to  modulo  2n  (integral  of  instantaneous  fre¬ 
quency  to  modulo  1),  and  since  the  modulo  of  a  sum  of  terms  is  equal  to  the  modulo  of  the  sum  of 
the  modulos  of  each  term,  let 


Zn  =  Z'n  modulo  1  =  Z„_!  +  — 

J  S 


modulo  1 . 


Applying  the  recurrence  formula  to  the  harmonically  related  instantaneous  frequencies  given  in 
Eq.  (lb)  and  using  the  parameters  for  the  <7 -subscripted  (or  higher  frequency)  signal,  gives 

Zq  n  =  Frac(Z?  „  _j  +  (fq  +  vn)/fs )  (2a) 


Zp,„  =  Frac{Zpn_,  +  (fq  +  vn)/qfs).  (2b) 

Here  Zp  =  Zq  =0,  and  the  function  FracJ. . .]  implies  taking  only  the  fractional  part  of  the 
decimal  argument.  It  is  recognized  that  /  =  qfp  is  the  mean  frequency  of  the  higher  frequency  sig¬ 
nal  and  vn  =  vq{nAt)  -  qvp{n  At)  is  the  fluctuating  signal  component  of  the  higher  frequency  sig¬ 
nal.  The  resulting  modulo-phase  functions  of  the  two  harmonically  related  signals  are 

ip  n  =  2irZ„  „  and  \pqn  =  2irZq  „  .  (2c) 

FORMULATION  OF  THE  FREQUENCY-FLUCTUATION  FUNCTION 

At  this  point,  a  particular  zero-mean  deterministic  function  for  may  be  chosen.  In  this 
report,  however,  the  frequency  function  is  chosen  to  be  truly  random  and  Gaussian-like  but  bounded 
in  its  peak  excursion  from  zero.  To  accomplish  this,  let  be  a  zero-mean  Gaussian  statistic  with 
standard  deviation  a ^  Hz.  To  bound  the  frequency-fluctuation  function  to  absolute  values  no  greater 
than  a  a r,  let 

vn  =  a  a-  Frac  ff„/a  or| ,  (3a) 

where  Frac)  |  implies  the  fractional  part  of  the  argument  while  retaining  the  sign  of  .  The  statis¬ 
tic  vn  is  therefore  a  zero-mean  random  number  whose  probability  density  is 

foo  .  ,  1 

2J  Pit  I  "I  +jaaf)  "  n'^ 

pv(»)  =  <  j  =0  -  (3b) 

0  for  o  <  |  v  | 


V 


$ 
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and 


P{(x)  = 


2a; 


V5*f 


The  probability  that  —  k0  <  p  <  p0  for  0  <  v0  <  a  is 


oo 


P(  |  V  \  <V0)  =  2  L  +  j<*  I  -  N[ja\], 

7=0 


where  N[y\  is  the  normal  or  Gaussian  distribution  function 

=  2  dx- 


(3c) 


FORMULATION  OF  THE  GAUSSIAN  STATISTIC  f„ 

To  emulate  the  frequency  fluctuations  of  real  signals,  the  rate  of  variation  (autocorrelation  or 
power  spectral  density)  of  the  statistic  must  be  controlled.  This  may  be  accomplished  by  low-pass 
filtering  a  running  sequence  of  independent  random  samples.  Consider,  then,  that  £„  is  a  zero-mean 
sequence  of  independent  random  samples  with  standard  deviation  af.  The  desired  output  can  be 
obtained  by  processing  £„  through  the  discrete  equivalence  of  a  single-pole  infinite-impulse  response 
(HR)  filter  [1,21,  whose  normalized  unit-impulse  response  is 

h„  =  (1  -  0)0"/,  ,  (4) 


where  /,  =  1/At  is  the  sample  rate  and  0  <  1  is  the  filter  design  parameter  that  controls  the  effec¬ 
tive  smoothing  time  or  noise  bandwidth  of  the  filter.  The  effective  bandwidth  of  the  discrete  filter  is 


B 


1  00 

7  E  hn  Al 

1  n=0 


1  ~  0  f_J_ 
1+02' 


(5a) 


and  the  effective  smoothing  time  is 


_L  =  1  +  <3  J_ 

2  B  1  ~  <3  f/ 


(5b) 


The  normalized  filter  autocorrelation  function  is,  from  Eqs.  (4)  and  (5b), 

In  =  T  £  h,h,i  ]n{At  =  =  0fA"  (50 

In  terms  of  the  smoothing  time  r,  the  bandwidth  B,  or  the  autocorrelation  at  a  time  delay 
t  (n  =  I  /At  =  /,/),  the  filter  parameter  becomes 
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„  _  fsT  ~  1  _  A-2  6  _  u-f, 

ft  r  *  ^  Tn  (6a) 

/v  r  +  1  fs  +  25 

The  analog  equivalent  of  this  discrete  filter  is  a  single-pole  low-pass  filter  whose  cutoff  frequency  (  3 
dB  down  point)  is 

,  1  25  1-5  fs 

/,=•—  =  —  =  7—7  —  •  (6b) 

XT  X  1  +  p  X 

The  effective  bandwidth  of  the  random  statistic  £„ ,  sampled  at  a  rate  /, ,  is  the  Nyquist  rate 
/,/ 2.  Consequently,  the  variance  of  £„  at  the  output  of  the  filter  is  reduced  by  the  factor 
2 B/fs  =  WfsT.  To  make  the  standard  deviation  at  the  filter  output  equal  to  the  desired  standard 
deviation  af,  the  random  statistic  £„  must  be  multiplied  by  the  factor  c  where,  from  Eq.  (5b), 


C  =  V/7  r  ^ 
ai 


l  +  5  £r 
l  -  ft 


The  output  of  the  filter  is  obtained  by  using  the  discrete  convolution  relation  [1,21 

?n  =  E  hJCZn  -]&>  ■ 

J  =0 

From  Eqs.  (4),  (6c)  and  (7a),  the  Gaussian  statistic  is  formulated  as 

=  (i  -  ft)  =  5rn-,  +  (i  -  ft)ctn 


ft? +  ^1"  -  52  — 


The  random  sequence  need  not  be  Gaussian  for  to  be  Gaussian  as  long  as 
/,  t  =  (I  +  5)/G  -  ft)  is  large  compared  with  one  (by  reason  of  the  central-limit  theorem).  The 
value  fsT  is  the  effective  number  of  independent  samples  of  £„  in  the  filter  smoothing  process.  This 
value  is  expected  to  be  very  large  in  practical  applications. 

FORMULATION  OF  THE  RANDOM  SEQUENCE 

A  convenient  method  of  generating  the  random  sequence  £„  on  the  digital  computer  is  to  ran 
domly  select  digital  values  Rr  between  0  and  1  at  the  sample  rate  J\  and  to  formulate  £„  as 

=  5„  -  0.5.  (Sat 

where  it  is  assumed  that  all  values  of  5„  are  equally  likely  to  occur  with  each  sample  and  thus  are 
independent.  The  random  statistic  has  zero  mean,  and  foi  M  discrete  values  between  0  and  1  the 
standard  deviation  of  is 


m 
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M  +  1 
12  (M  -  1) 


For  M  large,  the  standard  deviation  may  be  approximated  as  1  /  \AT2. 


SIMULATOR  DESIGN  EQUATIONS 

The  design  relations  to  simulate  two  narrowband  random  signals  embedded  in  random  noise  is 
now  given  in  an  orderly  manner  for  convenient  use  in  applications.  The  combined  signal  is  con¬ 
sidered  to  be  a  discrete  signal  sampled  at  a  rate  fs  and  expressed  in  the  form 

sn  =  ap  sin  (2t rZp  n)  +  aq  sin  (2irZq  „)  -I-  ,  (9) 

where  r\n  is  an  independent  zero-mean  Gaussian  statistic  with  variance  o and  the  functions  lie Zp  n 
and  2trZq  n  are  harmonically  related  random  phase  functions  (modulo  2ir). 


Input  Parameters 

The  following  parameters  are  to  be  entered  into  the  simulator  program. 


fs 

q>l 

Iq  =  lip 


P  =  r„  -  r 


is  the  signal  sample  rate  in  Hz, 

is  the  harmonic  ratio  of  the  narrowband  components, 

is  the  mean  frequency  of  the  higher  frequency  component  in  Hz 

is  the  standard  deviation  of  the  Gaussian  statistic  in  Hz 

is  the  standard  deviation  of  the  random  statistic  £„  in  Hz, 

is  the  standard  deviation  of  the  random-noise  statistic  T]n  , 

is  ratio  of  the  peak-frequency  fluctuation  of  vn  to  af, 

is  frequency-fluctuation  smoothing  time  in  seconds, 

is  the  <7 -channel  signal-to-noise  rms  ratio  in  dB/Hz,  and 

is  the  difference  in  the  p -channel,  and  q -channel  signal-to-noise 

rms  ratios  in  dB/Hz. 


Random  Signal  Generation 

The  first  step  in  signal  simulation  is  to  generate  two  random  sequences  of  independent  samples 
r)n  and  The  statistic  jj„  is  a  zero-mean  Gaussian  statistic  with  a  standard  deviation  of  av.  This 
statistic  is  used  as  the  broadband  noise  signal  in  Eq.  (9),  and  its  power  spectral  density  is  2 o~/J\. 
The  statistic  £„  is  a  zero-mean  random  statistic  with  a  standard  deviation  of  at.  This  statistic  is  used 
as  the  kernal  signal  in  the  generation  of  the  random  phase  functions  for  the  harmonically  related  com¬ 
ponents  in  Eq.  (9).  The  method  of  generating  the  two  random  sequences  are  left  to  the  user. 


%  **»  _  »  »  i'*  _»  _  _  »  »  «  ■  ij.%  a  «  .  d'*  o'*  jl^  o'*  k  "  k  *  *  *  • 


V  A 


GERLACH.  KUNZ.  ANDERSON.  AND  FLOWERS 


Program  Algorithms 

The  following  algorithm  formulations  simulate  the  discrete  signal  given  in  Eq.  (9): 


2  mV20 

a'  =  W'°m 


(10) 


ap  =  aq  10p/2°  =  -jjr  a,10' 


rj  20 


(ID 


0  = 


fs*  ~  1 

Is*  +  1 


(12) 


r„  =  /3r„_i  +  vi  -  02  —  ^ 


(13) 


for  n  =  0,  1 ,  2,  3,  .  .  . 


vn  =  Q-0f  Fracjfn/o;  af|  (14) 

where  Frac(-  ■  |  implies  the  decimal  fractional  part  of  the  argument  while  retaining  the  sign  of  the 
argument. 


=  Frac(Z9„_| 

+  if q 

+  vn)/fs\ 

(Zq ,  - 1 

=  0) 

(15) 

=  Frac  [Zpn_] 

+  (j fq 

+  vnVqfs  1 

(Zp.-i 

=  0) 

(16) 

sn  =  ap  sin  (2irZp  „)  +  aq  sin  (2irZq  n)  +  Vn  .  (17) 

EXAMPLES 

To  exemplify  the  use  of  the  simulator,  the  following  design  parameters  are  selected  to  demon¬ 
strate  the  performance  characteristics. 

/,  =  256  Hz  sample  rate, 

q  =1.75  harmonic  ratio, 

fq  =  21  Hz  </“  =  12  Hz)  , 
o (■  =  0.125  Hz. 

a  =3,  and 

r  =  36  s  (effective  smoothing  time). 

These  parameters  imply  that  the  upper  signal  frequency  variation  is  hounded  between  21  ±  0.375 
Hz,  and  the  e -folding  time  of  the  fluctuating-frequency  component  is  18s.  (The  e -folding  time  is  the 
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time  that  the  autocorrelation  of  the  fluetuating-frequency  component  is  equal  to  1/e.)  For  harmoni¬ 
cally  related  signals,  the  instantaneous  frequency  of  the  upper  signal  is  precisely  q  times  the  instan¬ 
taneous  frequency  of  the  lower  signal. 

Spectrogram  of  Example  Signal  Simulations 

Harmonically  related  signals  using  the  above  parameters  are  generated  along  with  a  random 
noise  signal.  The  signal-to-noise  spectral  density  ratio  of  both  signal  components  was  set  to  8  dB/Hz. 
A  spectrogram  of  the  resulting  signal  is  shown  in  Fig.  1.  The  resolution  of  the  spectrogram  in  the 
example  is  1/16  Hz.  The  two  harmonically  related  signal  components  are  located  at  center  frequen¬ 
cies  of  12  and  21  Hz  respectively.  Superimposed  on  the  spectrogram  for  comparison  is  a  plot  of  the 
actual  fluctuating  frequency  v„  given  by  Eq.  (14).  The  scale  of  this  plot  was  adjusted  to  closely 
approximate  that  of  the  upper-frequency  signal  and  placed  over  the  spectrogram  plot. 


10  20  30 

FREQUENCY  (Hz) 

Fig.  1  -  Spectrogram  of  simulated  harmonic  signals  in  broadband  noise.  The  mean  frequencies  of  the 
signals  are  12  and  21  Hz  (q  =  1.75),  and  the  signal-to-noise  ratio  of  each  signal  is  8  dB/Hz.  The  reso¬ 
lution  of  the  spectrogram  is  1/16  Hz  Superimposed  on  the  spectrogram  for  comparison  lat  approxi¬ 
mately  26  Hz)  is  a  plot  of  the  actual  fluctuating  frequency  components  v.  Parameters  of  the  signal  are 
given  in  the  accompanying  text. 

Phase  Characteristics  of  the  Simulated  Signals 

The  phase  characteristics  of  the  harmonically  related  signals  are  obtained  by  using  Hanning- 
windowed  sectionalized  Fourier  transforms  (SFTs)  to  Filter,  baseband,  and  decimate  the  two  nar¬ 
rowband  digital  signals.  The  SFT  sizes  (integration  times)  for  the  upper  and  lower  frequency  signals 
are  2  s  and  3.5  s,  respectively,  to  contain  the  spectrum  of  the  signals  within  the  main  filter  passbands 
of  0.75  Hz  and  3/7  Hz  respectively  [3|.  The  signal-to-noise  power  spectral  density  ratio  in  this  case 
is  30  dB  /Hz.  The  amplitude  and  phase  characteristics  at  the  filtered  outputs  over  a  10-min.  interval 
are  shown  in  Fig.  2. 

The  top  two  graphs  (Fig.  2)  are  the  plots  of  the  amplitude  and  phase  (harmonically  transformed) 
of  the  lower-frequency  signal,  and  the  middle  two  graphs  are  similar  plots  for  the  upper-frequency 
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0  123456789  10 

TIME  (min) 

Pig  2  Amplitude  and  phase  plots  of  simulated  harmonic  signals  showing  the  coherent 
nature  of  harmonically  related  signals  The  signal  to  noise  ratio  of  the  two  signals  is  30 
dB/H/,  and  the  signals  have  been  processed  through  proportional  SPTs  whose  bandwidths 
encompass  the  spectrum  of  the  signals  The  lop  two  graphs  are  temporal  plots  of  the 
amplitude  and  phase  of  the  lower-frequency  signal,  and  the  middle  two  graphs  are  similar 
plots  of  the  upper  frequency  signal  The  bottom  two  graphs  show  the  amplitude  product  and 
the  phase  difference  of  the  two  signals.  Harmonic  correlation  of  the  signals  is  the  integration 
of  the  amplitude  product  times  the  cosine  of  the  phase  difference  over  a  selected  time 
interval 

signal.  The  observed  amplitude  fluctuations  are  the  result  of  the  nonuniform  filter  characteristics 
(scallop  effect)  of  the  SFT  process  as  the  instantaneous  frequencies  vary  over  time.  The  similarity  of 
the  harmonically  related  signals  is  quite  evident.  The  bottom  two  graphs  are  plots  of  the  product  of 
the  two  amplitudes  and  the  phase  difference  between  the  two  signals.  This  is  of  interest  since  the  har¬ 
monic  correlation  between  the  two  signals  is  the  integration  of  the  amplitude  product  times  the  cosine 
of  the  phase  difference  over  a  given  smoothing  time.  Because  the  phase  difference  is  essentially  zero, 
the  harmonic  correlation  is  effectively  equal  to  the  integral  of  the  priniuct  of  the  amplitudes  as 
expected  for  harmonically  related  narrowband  signals  (The  noise  in  me  two  separated  filter  bands  is 
for  all  practical  purposes  uncorrelated) 
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To  demonstrate  the  influence  of  the  fiuctuating-frequency  components  vn  on  harmonic  correla¬ 
tion,  the  above  experiment  was  repeated  by  using  independent  and  uncorrelated  fiuctuating-frequency 
components  for  the  two  narrowband  signals.  The  mean  frequencies  and  the  standard  deviations  of  the 
fluctuating-  signal  components  remain  unchanged.  Figure  3  shows  the  results. 

The  format  of  the  plots  in  Fig.  3  is  the  same  as  for  Fig.  3.  In  this  case,  however,  there  is  little 
similarity  between  the  amplitudes  and  phases  of  she  two  signals.  In  particular,  the  phase  difference 
between  the  two  signals  (bottom  graph)  is  radically  variable  over  time.  Thus,  over  a  reasonable  time 
increment,  the  integral  of  the  amplitude  product  and  the  cosine  of  the  phase  difference  can  be 
expected  to  be  near  zero,  resulting  in  a  low  (or  zero)  correlation.  This  is  to  be  expected  for  nar¬ 
rowband  signals  that  are  not  harmonically  related. 
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Fig.  3  -  Amplitude  and  phase  plots  of  simulated  signals  illustrating  the  incoherent  nature  of 
the  signals  when  their  fluctuating  frequency  components  are  independent  and  uncorrelated. 
The  format  of  the  illustration  and  the  signal  parameters  are  the  same  as  in  Fig.  2  with  the  one 
exception  that  their  random  fluctuating  frequency  components  are  independent.  The  phase 
difference  between  the  signals  is  seen  to  vary  radically  over  time,  resulting  in  a  low  (or  zero) 
correlation  over  relatively  long  integration  intervals  (greater  than  I  min). 
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Instantaneous  Frequency  Distribution 

The  fluctuating  frequency  v  and  its  distribution  (or  probability  density)  are  given  in  Eq.  13),  To 
verify  the  distribution,  approximately  4  million  independent  samples  of  v  were  accumulated,  and  the 
distribution  of  via f  was  computed.  The  results  are  shown  in  Fig.  4. 

The  data  in  Fig.  4  are  plotted  as  a  function  of  the  normalized  variable  via f  for  a  =  3.  The 
continuous  curve  is  the  theoretical  probability  density  given  by  Eq.  (3b).  The  dots  show  the  results 
computed  from  the  experimental  data.  Although  these  results  (computed  from  the  large  sample  size) 
closely  approximate  the  theoretical  statistical  distribution,  it  can  be  expected  that  results  for  a  small 
sample  size  (such  as  that  realized  over  a  period  of  30  min)  will  deviate  significantly  from  the  theoreti¬ 
cal  predictions. 


Via* 

Fig  4  —  Normalized  instantaneous  frequency  distribution  of  the  simulated 
fluctuating-frequency  component  for  a  =  3.  The  continuous  curve  is  the 
theoretical  probability  density,  and  the  solid  dots  are  the  experimental  results 
computed  from  approximately  4-million  independent  samples  of  the  fluctuat¬ 
ing  frequency  function  v.  For  a  small  sample  size,  the  experimental  distri¬ 
bution  can  deviate  significantly  from  the  theoretical  curve. 

CONCLUSIONS 


1 .  A  relatively  simple  algorithm  is  formulated  to  simulate  harmonically  related  narrowband  signals 
with  random  frequency  fluctuations  in  discrete  format  on  a  digital  computer. 

2.  Signal  parameters  are  incorporated  into  the  algorithm  to  select  the  signal  mean  frequencies  and 
to  control  both  the  spectral  bounds  and  the  autocorrelation  (or  power  spectral  density)  of  the  sig¬ 
nal  fluctuations. 

3.  Examples,  using  the  algorithm,  demonstrate  the  performance  of  the  signal  simulator  and  its  con¬ 
formance  with  theoretical  predictions.  The  results  are  shown  in  Figs.  I  through  4. 
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